function [ correctedImageSet  scoreMatrix  magUniformCorrection bestBlueFrac bestOtherFrac ] = ...
    computeUniformBgCorrection(imageSet, imageFeatures, loPercentile)

% % Below code is for when we were not optimizing the uniform background
% % correction. It utilizes the runUniformBackgroundCorrection function
% loPercentileFraction = 1;
% [correctedImageSet.blueIm, magBlueCorrection] = runUniformBackgroundCorrection(imageSet.blueIm,...
%     loPercentile,loPercentileFraction);
% [correctedImageSet.otherIm, magOtherCorrection] = runUniformBackgroundCorrection(imageSet.otherIm,...
%     loPercentile,loPercentileFraction);
% 
% % Output the original red image in the redIm field for struct.field consistancy
% correctedImageSet.redIm = imageSet.redIm;
% 
% magUniformCorrection.blueCenter = magBlueCorrection;
% magUniformCorrection.blueBoundary = magBlueCorrection;
% magUniformCorrection.otherCenter = magOtherCorrection;
% magUniformCorrection.otherBoundary = magOtherCorrection;
% 
% bestScore = -1;

% 

% First determine the lowest 1/2 percentile pixel values. This is our
% original uniform background correction, and we want to optimize by taking
% fractions

% Compute nonZeroImVals for both blue and other
blueNonZeroImVals = imageSet.blueIm(imageSet.blueIm ~= 0);
otherNonZeroImVals = imageSet.otherIm(imageSet.otherIm ~= 0);

% Define the hiPercentile based on the loPercentile input
hiPercentile = 100 - loPercentile;

blueLoPercentileVal = prctile(blueNonZeroImVals, loPercentile);
otherLoPercentileVal = prctile(otherNonZeroImVals, loPercentile);

% Now generate a list of fractions to multiply the loPercentile value.
% Perform a 2-D grid search for best fraction

blueRangeFractionVals = 0:0.1:1;
otherRangeFractionVals = 0:0.1:1;

% Flip dims so we start at the original low value
blueFractionVals = flipdim(blueRangeFractionVals, 2);
otherFractionVals = flipdim(otherRangeFractionVals, 2);

% Write 2-D grid search similar to variable correction
% computeRedBgCorrection function. We want to use the same processing
% techniques to allow for the most accurate comparision between methods.
bestBlueFrac = 1;
bestOtherFrac = 1;
bestRatioScore = 10^6; % initialize to worst possible variance value
scoreMatrix = zeros(length(blueFractionVals), length(otherFractionVals));
i = 1;
j = 1;
rowRange = imageFeatures.fovRow - imageFeatures.macRadius: imageFeatures.fovRow + imageFeatures.macRadius;
colRange = imageFeatures.fovCol - imageFeatures.macRadius: imageFeatures.fovCol + imageFeatures.macRadius;

for blueFrac = blueFractionVals
    for otherFrac = otherFractionVals
        blueCorrected = imageSet.blueIm - blueFrac*blueLoPercentileVal;
        otherCorrected = imageSet.otherIm - otherFrac*otherLoPercentileVal;
        blueCorrected(imageFeatures.iNonMacRegion) = mean(blueCorrected(imageFeatures.iMacBoundary));
        otherCorrected(imageFeatures.iNonMacRegion) = mean(otherCorrected(imageFeatures.iMacBoundary));
        ratioIm = blueCorrected ./ otherCorrected;
        ratioIm(imageFeatures.iNonMacRegion) = ...
            mean(ratioIm(imageFeatures.iMacBoundary));
        % Ignore ratios with more than 10 negative pixel values. These are
        % deemed implausible results
        negValueLocations = any(ratioIm < 0);
        numNegValues = sum(negValueLocations);
        if numNegValues > 10
            ratioScore = 100000;
        else
            ratioScore = compute_total_variation(ratioIm(rowRange, colRange));
        end
        
        scoreMatrix(i,j) = ratioScore;
        j = j+1;
        if ratioScore < bestRatioScore
            bestRatioScore = ratioScore;
            bestBlueFrac = blueFrac;
            bestOtherFrac = otherFrac;
        end
    end
        i = i+1;
        j = 1;
end

% Compute best corrected image and renormalize
correctedImageSet.blueIm = (imageSet.blueIm - bestBlueFrac*blueLoPercentileVal)...
    /(1-2*bestBlueFrac*blueLoPercentileVal);
correctedImageSet.otherIm = (imageSet.otherIm - bestOtherFrac*otherLoPercentileVal)...
    /(1-2*otherFrac*otherLoPercentileVal);
correctedImageSet.redIm = imageSet.redIm;

magUniformCorrection.blueCenter = bestBlueFrac*blueLoPercentileVal;
magUniformCorrection.otherCenter = bestOtherFrac*otherLoPercentileVal;  
magUniformCorrection.blueBoundary = magUniformCorrection.blueCenter;
magUniformCorrection.otherBoundary = magUniformCorrection.otherCenter;

end

